% this file is to be run after estimation is over
% and all the estimated files are downloaded from server
% this includes ESTIMATES.mat in this folder and 
% theta_star_estimates in /Bootstrap_CI/BS_theta_estimates
clear
clc
close all 
thisFile = mfilename('fullpath');
thisDir = fileparts(thisFile);
cd(fileparts(fileparts(thisDir)));
%% Table 3
load('output/estimates/ESTIMATES.mat')
% order 
% theta_mle      = [param_est(1), ... % beta
%                   param_est(2), ... % lambda
%                   param_est(3), ... % muU
%                   param_est(4), ... % sig^2U
%                   param_est(5:ng+4),... %alv
%                   param_est(ng+5:2*ng+4), ... % betv
%                   param_est(2*ng+5), ... %a1-tilde
%                   param_est(2*ng+6), ... %a2-tilde
%                   param_est(2*ng+7), ... %bar-w
%                   param_est(end), ... %ulow
%                   ];
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% confidence interval 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    BS_size =150; % number of bootstrap sample
    BS_sample=zeros(BS_size,size(theta_mle,2)) ;
    filename = 'output/bootstrap/BS_theta_estimates/theta_star_estimates_';
    for kk=1:BS_size
        Filename=strcat(filename, num2str(kk), '.mat');
        BS_sample(kk,:)=struct2array(load(Filename));
    end
    CI = zeros(size(theta_mle,2),2); %[2.5% and 97.5%]
    for kk =1:length(CI)
        CI(kk,:)=prctile(BS_sample(:,kk),[2.5,97.5]);
    end

    CI = array2table(CI, "VariableNames",{'2.5%', '97.5%'});

save('output/tables_and_figures/CI.mat', 'CI')
%% Figure 5 Estimated PDFs of Costs and Demand Shocks
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfig1: U~f_U = Truncated Normal (mu_u, sigma_u^2) truncated at u_low
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f_u = makedist('Normal','mu',theta_mle(3),...
    'sigma',sqrt(theta_mle(4))); % we estimate variance but matlab takes std. dev. 
f_u = truncate(f_u, theta_mle(20),Inf);
U = linspace(floor(theta_mle(end)/10)*10,ceil((theta_mle(3)...
    + 2.5*sqrt(theta_mle(4)))/10)*10,1000);

subplot(1,3,1)
plot(U,pdf(f_u,U), "LineWidth",3)
grid on 
title('f_{U^{dt}}(\cdot)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%subfig2: f_{W|U}(.| U = {25%, 50%, 75%}percentile)
% W|U=u ~ w_bar x (2x Beta(a_w(u), a_w(u))-1) where 
% a_w(u) = exp(a_1_tilde + a_2_tilde x u) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find 25th, 50th, 75th percentile of U for conditioning
u_draws = random(f_u,100000,1);
u_percentiles = prctile(u_draws,[25, 50,75]); 


a_w_u = @(u) exp(theta_mle(17)+theta_mle(18)*u);

% define the densities 
f_w_u_25=@(w)betapdf((w+theta_mle(19))/theta_mle(19)/2,a_w_u(u_percentiles(1))...
    ,a_w_u(u_percentiles(1)))/2/theta_mle(19); 
f_w_u_50=@(w)betapdf((w+theta_mle(19))/theta_mle(19)/2,a_w_u(u_percentiles(2))...
    ,a_w_u(u_percentiles(2)))/2/theta_mle(19);
f_w_u_75=@(w)betapdf((w+theta_mle(19))/theta_mle(19)/2,a_w_u(u_percentiles(3))...
    ,a_w_u(u_percentiles(3)))/2/theta_mle(19);



%plot the densities 
subplot(1,3,2)
hold on 
plot(linspace(-25,25,1000),f_w_u_25(linspace(-25,25,1000)), 'LineWidth',3)
plot(linspace(-25,25,1000),f_w_u_50(linspace(-25,25,1000)), 'LineWidth',3, 'LineStyle','--')
plot(linspace(-25,25,1000),f_w_u_75(linspace(-25,25,1000)), 'LineWidth',3, 'LineStyle','-.')
hold off 
grid on 
title('f_{W|U^{dt}}(\cdot|\cdot)')
legend({'U=p25', 'U=p50', 'U=p75'}, "Location","best")


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%subfig3: V~f_{V}= w_bar x (Beta(a_i, b_i)+1)
% where i=1,...6 groups. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% define the density parameters
ng =6; % number of groups
f_v = cell(ng,1); %define density 
a= theta_mle(5:ng+4); % Beta parameter a
b = theta_mle(ng+5:2*ng+4);% Beta parameter b
v_space = linspace(theta_mle(19),2*theta_mle(19),5000);

for kk=1:6
    f_v{kk} = @(v)betapdf(v./theta_mle(19)-1,a(kk),b(kk))/theta_mle(19); %define pdf f_V
end

% plot the densities 
subplot(1,3,3)
hold on
for kk=1:ng
    plot(v_space,f_v{kk}(v_space), 'LineWidth',2)
end
hold off
grid on
title('f_{V}(\cdot)')
legend({'Group 1', 'Group 2', 'Group 3','Group 4','Group 5', 'Group 6'}, "Location","best")


saveas(gcf,'output/tables_and_figures/pdfshocks','epsc')
saveas(gcf,'output/tables_and_figures/pdfshocks','fig')


